Tutorial: Running populations with binary_c-python (grid sampling)
This notebook will show you how to evolve a population of stars
Much of the code in the binary_c-python package is written to evolve a population of stars through the Population object, rather than evolving a single system. Let’s go through the functionality of this object step by step and set up some example populations.
We will use a grid-based sampling technique in this notebook. For Monte-Carlo based sampling, see this notebook
At the bottom of this notebook there are some complete example scripts.
[1]:
import os
from binarycpython.utils.functions import temp_dir, output_lines
from binarycpython import Population
TMP_DIR = temp_dir("notebooks", "notebook_population", clean_path=True)
# help(Population) # Uncomment to see the public functions of this object
Setting up the Population object
To set up and configure the population object we need to make an object instance of the Population object, and add configuration via the .set() function.
There are three categories of options that the Population object can set:
BSE options: these options will be used for the binary_c calls, and are recognized by comparing the arguments to a known list of available arguments of binary_c. To see which options are available, see section
binary_c parametersin the documentation. You can access these throughpopulation.bse_options['<bse option name>']after you have set them.Grid options: these options will be used to configure the behaviour of the Population object. To see which options are available, see section
Population grid code optionsin the documentation. They can be accessed viapopulation.population_options['<grid option name>']after you have set them.Custom options: these options are not recognized as either of the above, so they will be stored in the custom_options, and can be accessed via
population.custom_options['<custom option name>']
[2]:
# Create population object
example_pop = Population(tmp_dir=TMP_DIR)
# If you want verbosity, set this before other things
example_pop.set(verbosity=1)
# Setting values can be done via .set(<parameter_name>=<value>)
# Values that are known to be binary_c_parameters are loaded into bse_options.
# Those that are present in the default population_options are set in population_options
# All other values that you set are put in a custom_options dict
example_pop.set(
# binary_c physics options
M_1=10, # bse_options
orbital_period=45000000080, # bse_options
max_evolution_time=15000, # bse_options
eccentricity=0.02, # bse_options
# population_options
num_cores=2,
# Custom options # TODO: need to be set in population_options probably
data_dir=os.path.join(
TMP_DIR, "example_python_population_result"
), # custom_options
base_filename="example_pop.dat", # custom_options
)
# We can use the options through
print(example_pop.population_options['verbosity'])
print(example_pop.custom_options['base_filename'])
print(example_pop.bse_options['M_1'])
1
example_pop.dat
10
After configuring the population, but before running the actual population, its usually a good idea to export the full configuration (including version info of binary_c and all the parameters) to a file. To do this we use example_pop.export_all_info().
On default this exports everything, each of the sections can be disabled:
population settings (bse_options, population_options, custom_options), turn off with include_population settings=False
binary_c_defaults (all the commandline arguments that binary c accepts, and their defaults). turn off with include_binary_c_defaults=False
include_binary_c_version_info (all the compilation info, and information about the compiled parameters), turn off with include_binary_c_version_info=False
include_binary_c_help_all (all the help information for all the binary_c parameters), turn off with include_binary_c_help_all=Fase
On default it will write this to the custom_options[‘data_dir’], but that can be overriden by setting use_datadir=False and providing an outfile=<>
[3]:
example_pop.export_all_info()
[3]:
'/tmp/binary_c_python-david/notebooks/notebook_population/example_python_population_result/example_pop_settings.json.gz'
Adding grid variables
The main purpose of the Population object is to handle the population synthesis side of running a set of stars. The main method to do this with binary_c-python, as is the case with Perl binarygrid, is to use grid variables. These are loops over a predefined range of values, where a probability will be assigned to the systems based on the chosen probability distributions.
Usually we use either 1 mass grid variable, or a trio of mass, mass ratio and period (See below for full examples of all of these). We can, however, also add grid sampling for e.g. eccentricity, metallicity or other parameters.
In some cases it could be easier to set up a for loop that sets that parameter and calls the evolve function several times, e.g. when you want to vary a prescription (usually a discrete, unweighted parameter)
A notable special type of grid variable is that of the Moe & di Stefano 2017 dataset (see further down in the notebook).
To add a grid variable to the population object we use example_pop.add_sampling_variable (see next cell)
[4]:
help(example_pop.add_sampling_variable)
Help on method add_sampling_variable in module binarycpython.utils.population_extensions.sampling_variables:
add_sampling_variable(name: str, parameter_name: str, longname: str, valuerange: Union[list, str], samplerfunc: str, probdist: str, dphasevol: Union[str, int] = -1, gridtype: str = 'centred', branchpoint: int = 0, branchcode: Optional[str] = None, precode: Optional[str] = None, postcode: Optional[str] = None, topcode: Optional[str] = None, bottomcode: Optional[str] = None, condition: Optional[str] = None, dry_parallel: Optional[bool] = False, dependency_variables: Optional[list] = None) -> None method of binarycpython.utils.population_class.Population instance
Function to add sampling variables to the population_options.
The execution of the sampling generation will be through a nested for loop.
Each of the grid variables will get create a deeper for loop.
The real function that generates the numbers will get written to a new file in the TMP_DIR,
and then loaded imported and evaluated.
beware that if you insert some destructive piece of code, it will be executed anyway.
Use at own risk.
Args:
name:
name of parameter used in the grid Python code.
This is evaluated as a parameter and you can use it throughout
the rest of the function
Examples::
name = 'lnM_1'
parameter_name:
name of the parameter in binary_c
This name must correspond to a Python variable of the same name,
which is automatic if parameter_name == name.
Note: if parameter_name != name, you must set a
variable in "precode" or "postcode" to define a Python variable
called parameter_name
longname:
Long name of parameter
Examples::
longname = 'Primary mass'
range:
Range of values to take. Does not get used really, the samplerfunc is used to
get the values from
Examples::
range = [math.log(m_min), math.log(m_max)]
samplerfunc:
Function returning a list or numpy array of samples spaced appropriately.
You can either use a real function, or a string representation of a function call.
Examples::
samplerfunc = "self.const_linear(math.log(m_min), math.log(m_max), {})".format(resolution['M_1'])
precode:
Extra room for some code. This code will be evaluated within the loop of the
sampling function (i.e. a value for lnM_1 is chosen already)
Examples::
precode = 'M_1=math.exp(lnM_1);'
postcode:
Code executed after the probability is calculated.
probdist:
Function determining the probability that gets assigned to the sampled parameter
Examples::
probdist = 'self.Kroupa2001(M_1)*M_1'
dphasevol:
part of the parameter space that the total probability is calculated with. Put to -1
if you want to ignore any dphasevol calculations and set the value to 1
Examples::
dphasevol = 'dlnM_1'
condition:
condition that has to be met in order for the grid generation to continue
Examples::
condition = "self.population_options['binary']==1"
gridtype:
Method on how the value range is sampled. Can be either 'edge' (steps starting at
the lower edge of the value range) or 'centred'
(steps starting at ``lower edge + 0.5 * stepsize``).
dry_parallel:
If True, try to parallelize this variable in dry runs.
topcode:
Code added at the very top of the block.
bottomcode:
Code added at the very bottom of the block.
dependency_variables:
TODO: describe
All the distribution functions that we can use are stored in the binarycpython.utils.distribution_functions or binarycpython/utils/distribution_functions.py on git. If you uncomment the help statement below you can see which functions are available now:
[5]:
# import binarycpython.utils.distribution_functions
# help(binarycpython.utils.distribution_functions)
The next cell contains an example of adding the mass grid variable, but sampling in log mass. The commented grid variables are examples of the mass ratio sampling and the period sampling.
[6]:
# Add grid variables
resolution = {"M_1": 20}
# Mass
example_pop.add_sampling_variable(
name="lnm1",
longname="Primary mass",
valuerange=[2, 150],
samplerfunc="self.const_linear(math.log(2), math.log(150), {})".format(resolution["M_1"]),
precode="M_1=math.exp(lnm1)",
probdist="self.three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
dphasevol="dlnm1",
parameter_name="M_1",
condition="", # Impose a condition on this grid variable. Mostly for a check for yourself
)
# # Mass ratio
# test_pop.add_sampling_variable(
# name="q",
# longname="Mass ratio",
# valuerange=["0.1/M_1", 1],
# samplerfunc="self.const_linear(0.1/M_1, 1, {})".format(resolution['q']),
# probdist="self.flatsections(q, [{'min': 0.1/M_1, 'max': 1.0, 'height': 1}])",
# dphasevol="dq",
# precode="M_2 = q * M_1",
# parameter_name="M_2",
# condition="", # Impose a condition on this grid variable. Mostly for a check for yourself
# )
# #
# test_pop.add_sampling_variable(
# name="log10per", # in days
# longname="log10(Orbital_Period)",
# valuerange=[0.15, 5.5],
# samplerfunc="self.const_linear(0.15, 5.5, {})".format(resolution["per"]),
# precode="""orbital_period = 10** log10per
# sep = calc_sep_from_period(M_1, M_2, orbital_period)
# sep_min = calc_sep_from_period(M_1, M_2, 10**0.15)
# sep_max = calc_sep_from_period(M_1, M_2, 10**5.5)""",
# probdist="self.sana12(M_1, M_2, sep, orbital_period, sep_min, sep_max, math.log10(10**0.15), math.log10(10**5.5), -0.55)",
# parameter_name="orbital_period",
# dphasevol="dlog10per",
# )
Added sampling variable: {
"name": "lnm1",
"parameter_name": "M_1",
"longname": "Primary mass",
"valuerange": [
2,
150
],
"samplerfunc": "self.const_linear(math.log(2), math.log(150), 20)",
"precode": "M_1=math.exp(lnm1)",
"postcode": null,
"probdist": "self.three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
"dphasevol": "dlnm1",
"condition": "",
"gridtype": "centred",
"branchpoint": 0,
"branchcode": null,
"topcode": null,
"bottomcode": null,
"sampling_variable_number": 0,
"dry_parallel": false,
"dependency_variables": null
}
Setting logging and handling the output
On default, binary_c will not output anything (except for ‘SINGLE STAR LIFETIME’). It is up to us to determine what will be printed. We can either do that by hardcoding the print statements into binary_c (see documentation binary_c). Or, we can use the custom logging functionality of binarycpython (see notebook notebook_custom_logging.ipynb), which is faster to set up and requires no recompilation of binary_c, but is somewhat more limited in its functionality.
After configuring what will be printed, we need to make a function to parse the output. This can be done by setting the parse_function parameter in the population object (see also notebook notebook_individual_systems.ipynb).
In the code below we will set up both the custom logging, and a parse function to handle that output
[7]:
# Create custom logging statement: in this case we will log when the star turns into a compact object, and then terminate the evolution.
custom_logging_code = """
if(stardata->star[0].stellar_type >= 13)
{
if (stardata->model.time < stardata->model.max_evolution_time)
{
Printf("EXAMPLE_COMPACT_OBJECT %30.12e %g %g %g %d\\n",
//
stardata->model.time, // 1
stardata->star[0].mass, // 2
stardata->common.zero_age.mass[0], // 3
stardata->model.probability, // 4
stardata->star[0].stellar_type // 5
);
};
/* Kill the simulation to save time */
stardata->model.max_evolution_time = stardata->model.time - stardata->model.dtm;
};
"""
example_pop.set(
C_logging_code=custom_logging_code
)
The parse function must now catch lines that start with “EXAMPLE_COMPACT_OBJECT”, and write that line to a file
[8]:
def parse_function(self, output):
"""
Example parse function
"""
# get info from the population instance
data_dir = self.custom_options["data_dir"]
base_filename = self.custom_options["base_filename"]
# Check directory, make if necessary
os.makedirs(data_dir, exist_ok=True)
seperator = " "
# Create filename
outfilename = os.path.join(data_dir, base_filename)
parameters = ["time", "mass", "zams_mass", "probability", "stellar_type"]
# Go over the output.
for line in output_lines(output):
headerline = line.split()[0]
# CHeck the header and act accordingly
if headerline == "EXAMPLE_COMPACT_OBJECT":
values = line.split()[1:]
print(line)
if not len(parameters) == len(values):
print("Number of column names isnt equal to number of columns")
raise ValueError
if not os.path.exists(outfilename):
with open(outfilename, "w") as f:
f.write(seperator.join(parameters) + "\n")
with open(outfilename, "a") as f:
f.write(seperator.join(values) + "\n")
# Add the parsing function
example_pop.set(
parse_function=parse_function,
)
Evolving the grid
Now that we configured all the main parts of the population object, we can actually run the population! Doing this is straightforward: example_pop.evolve()
This will start up the processing of all the systems. We can control how many cores are used by settings num_cores. By setting the verbosity of the population object to a higher value we can get a lot of verbose information about the run, but for now we will set it to 0.
There are many population_options that can lead to different behaviour of the evolution of the grid. Please do have a look at those: grid options docs, and try
[9]:
# change verbosity
example_pop.set(verbosity=0)
## Executing a population
## This uses the values generated by the grid_variables
analytics = example_pop.evolve() # TODO: update this function call
Grid has handled 19 stars with a total probability of 0.0443872
**********************************
* Dry run *
* Total starcount is 19 *
* Total probability is 0.0443872 *
**********************************
Signalling processes to stop
EXAMPLE_COMPACT_OBJECT 3.598268106262e+01 1.30592 8.75988 0.00193614 13
EXAMPLE_COMPACT_OBJECT 2.436983545180e+01 1.35842 10.9948 0.00144093 13
EXAMPLE_COMPACT_OBJECT 1.690157943855e+01 1.43124 13.7998 0.00107238 13
EXAMPLE_COMPACT_OBJECT 1.242397939734e+01 1.52416 17.3205 0.000798096 13
EXAMPLE_COMPACT_OBJECT 9.758088010684e+00 1.6691 21.7394 0.000593966 13
EXAMPLE_COMPACT_OBJECT 8.484238931453e+00 1.66351 27.2857 0.000442046 13
EXAMPLE_COMPACT_OBJECT 7.699271848690e+00 1.64817 34.247 0.000328983 13
EXAMPLE_COMPACT_OBJECT 7.568274926111e+00 1.64805 42.9844 0.000244839 13
EXAMPLE_COMPACT_OBJECT 7.569952912843e+00 1.64805 53.9508 0.000182216 13
EXAMPLE_COMPACT_OBJECT 7.571923249967e+00 1.64787 67.7151 0.00013561 13
EXAMPLE_COMPACT_OBJECT 7.612944004713e+00 1.64836 84.9909 0.000100925 13
EXAMPLE_COMPACT_OBJECT 7.620339349962e+00 1.64833 106.674 7.51114e-05 13
****************************************************
* Process 0 finished: *
* generator started at 2023-05-18T18:01:33.543282 *
* generator finished at 2023-05-18T18:01:34.555599 *
* total: 1.01s *
* of which 0.94s with binary_c *
* Ran 9 systems *
* with a total probability of 0.0234439 *
* This thread had 0 failing systems *
* with a total failed probability of 0 *
* Skipped a total of 0 zero-probability systems *
* *
****************************************************
EXAMPLE_COMPACT_OBJECT 7.622200536471e+00 1.64833 133.89 5.59e-05 13
****************************************************
* Process 1 finished: *
* generator started at 2023-05-18T18:01:33.547219 *
* generator finished at 2023-05-18T18:01:34.634907 *
* total: 1.09s *
* of which 1.02s with binary_c *
* Ran 10 systems *
* with a total probability of 0.0209432 *
* This thread had 0 failing systems *
* with a total failed probability of 0 *
* Skipped a total of 0 zero-probability systems *
* *
****************************************************
**********************************************************
* Population-365901e7ef5d46da9920af08b29842aa finished! *
* The total probability is 0.0443872. *
* It took a total of 1.53s to run 19 systems on 2 cores *
* = 3.05s of CPU time. *
* Maximum memory use 342.781 MB *
**********************************************************
No failed systems were found in this run.
After the run is complete, some technical report on the run is returned. I stored that in analytics. As we can see below, this dictionary is like a status report of the evolution. Useful for e.g. debugging.
[10]:
print(analytics)
{'population_id': '365901e7ef5d46da9920af08b29842aa', 'evolution_type': 'grid', 'failed_count': 0, 'failed_prob': 0, 'failed_systems_error_codes': [], 'errors_exceeded': False, 'errors_found': False, 'total_probability': 0.044387171445641534, 'total_count': 19, 'start_timestamp': 1684429293.5052047, 'end_timestamp': 1684429295.0318377, 'time_elapsed': 1.5266330242156982, 'total_mass_run': 649.905447944397, 'total_probability_weighted_mass_run': 0.28133908148630704, 'zero_prob_stars_skipped': 0}
Noteworthy functionality
Some extra features that are available from via the population object are:
write_binary_c_calls_to_file: Function to write the calls that would be passed to binary_c to a file
[11]:
help(example_pop.write_binary_c_calls_to_file)
Help on method write_binary_c_calls_to_file in module binarycpython.utils.population_extensions.dataIO:
write_binary_c_calls_to_file(output_dir: Optional[str] = None, output_filename: Optional[str] = None, include_defaults: bool = False, encoding='utf-8') -> None method of binarycpython.utils.population_class.Population instance
Function that loops over the grid code and writes the generated parameters to a file.
In the form of a command line call
Only useful when you have a variable grid as system_generator. MC wouldn't be that useful
Also, make sure that in this export there are the basic parameters
like m1,m2,sep, orb-per, ecc, probability etc.
TODO: this function can probably be cleaned a bit and can rely on the other startup and clean up functions (see population_class)
On default this will write to the datadir, if it exists
Args:
output_dir: (optional, default = None) directory where to write the file to. If custom_options['data_dir'] is present, then that one will be used first, and then the output_dir
output_filename: (optional, default = None) filename of the output. If not set it will be called "binary_c_calls.txt"
include_defaults: (optional, default = None) whether to include the defaults of binary_c in the lines that are written. Beware that this will result in very long lines, and it might be better to just export the binary_c defaults and keep them in a separate file.
Returns:
filename: filename that was used to write the calls to
[12]:
example_pop.set(verbosity=1)
calls_filename = example_pop.write_binary_c_calls_to_file()
print(calls_filename)
with open(calls_filename, 'r') as f:
print('\n'.join(f.read().splitlines()[:4]))
print("(abridged)")
Write grid code to /tmp/binary_c_python-david/notebooks/notebook_population/binary_c_grid_365901e7ef5d46da9920af08b29842aa.py [dry_run = False]
Grid has handled 19 stars with a total probability of 0.0443872
/tmp/binary_c_python-david/notebooks/notebook_population/example_python_population_result/binary_c_calls.txt
binary_c M_1 2.2406484012210224 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.22723621650191106 probability 0.011394572976608001
binary_c M_1 2.812296769855663 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.22723621650191117 probability 0.008480166685456411
binary_c M_1 3.5297876799548944 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.22723621650191106 probability 0.006311182276049824
binary_c M_1 4.430329401616038 eccentricity 0.02 max_evolution_time 15000 orbital_period 45000000080 phasevol 0.22723621650191106 probability 0.004696962123378559
(abridged)
Full examples of population scripts
Single star population
Below is a full setup for a population of single stars
[13]:
def parse_function(self, output):
"""
Example parsing function
"""
# extract info from the population instance
# Get some information from the
data_dir = self.custom_options["data_dir"]
base_filename = self.custom_options["base_filename"]
# Check directory, make if necessary
os.makedirs(data_dir, exist_ok=True)
#
seperator = " "
# Create filename
outfilename = os.path.join(data_dir, base_filename)
# The header columns matching this
parameters = ["time", "mass", "zams_mass", "probability", "radius", "stellar_type"]
# Go over the output.
for el in output_lines(output):
headerline = el.split()[0]
# CHeck the header and act accordingly
if headerline == "MY_STELLAR_DATA":
values = el.split()[1:]
if not len(parameters) == len(values):
print("Number of column names isnt equal to number of columns")
raise ValueError
if not os.path.exists(outfilename):
with open(outfilename, "w") as f:
f.write(seperator.join(parameters) + "\n")
with open(outfilename, "a") as f:
f.write(seperator.join(values) + "\n")
# Create population object
example_pop = Population()
# If you want verbosity, set this before other things
example_pop.set(verbosity=0)
# Setting values can be done via .set(<parameter_name>=<value>)
example_pop.set(
# binary_c physics options
M_1=10, # bse_options
separation=0, # bse_options
orbital_period=45000000080, # bse_options
max_evolution_time=15000, # bse_options
eccentricity=0.02, # bse_options
# population_options
num_cores=2,
tmp_dir=TMP_DIR,
# Custom options: the data directory and the output filename
data_dir=os.path.join(
TMP_DIR, "example_python_population_result"
), # custom_options
base_filename="example_pop.dat", # custom_options
)
# Creating a parsing function
example_pop.set(
parse_function=parse_function, # Setting the parse function thats used in the evolve_population
)
### Custom logging
# Log the moment when the star turns into neutron
example_pop.set(
C_logging_code="""
if(stardata->star[0].stellar_type >= 13)
{
if (stardata->model.time < stardata->model.max_evolution_time)
{
Printf("MY_STELLAR_DATA %30.12e %g %g %g %g %d\\n",
//
stardata->model.time, // 1
stardata->star[0].mass, // 2
stardata->common.zero_age.mass[0], // 4
stardata->model.probability, // 5
stardata->star[0].radius, // 6
stardata->star[0].stellar_type // 7
);
};
/* Kill the simulation to save time */
stardata->model.max_evolution_time = stardata->model.time - stardata->model.dtm;
};
"""
)
# Add grid variables
resolution = {"M_1": 20}
# Mass
example_pop.add_sampling_variable(
name="lnm1",
longname="Primary mass",
valuerange=[2, 150],
samplerfunc="self.const_linear(math.log(2), math.log(150), {})".format(resolution["M_1"]),
precode="M_1=math.exp(lnm1)",
probdist="self.three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
dphasevol="dlnm1",
parameter_name="M_1",
condition="",
)
# Exporting of all the settings can be done with .export_all_info()
example_pop.export_all_info()
# remove the result file if it exists
if os.path.isfile(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat")):
os.remove(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat"))
# Evolve the population
example_pop.evolve()
#
with open(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat"), 'r') as f:
output = f.read()
print("\n")
print(output)
Grid has handled 19 stars with a total probability of 0.0443872
**********************************
* Dry run *
* Total starcount is 19 *
* Total probability is 0.0443872 *
**********************************
Signalling processes to stop
****************************************************
* Process 0 finished: *
* generator started at 2023-05-18T18:01:37.291405 *
* generator finished at 2023-05-18T18:01:38.190228 *
* total: 0.90s *
* of which 0.83s with binary_c *
* Ran 9 systems *
* with a total probability of 0.0234439 *
* This thread had 0 failing systems *
* with a total failed probability of 0 *
* Skipped a total of 0 zero-probability systems *
* *
****************************************************
****************************************************
* Process 1 finished: *
* generator started at 2023-05-18T18:01:37.302546 *
* generator finished at 2023-05-18T18:01:38.250025 *
* total: 0.95s *
* of which 0.88s with binary_c *
* Ran 10 systems *
* with a total probability of 0.0209432 *
* This thread had 0 failing systems *
* with a total failed probability of 0 *
* Skipped a total of 0 zero-probability systems *
* *
****************************************************
**********************************************************
* Population-5b3001f6abf547c98152509881eadb3e finished! *
* The total probability is 0.0443872. *
* It took a total of 1.37s to run 19 systems on 2 cores *
* = 2.73s of CPU time. *
* Maximum memory use 350.297 MB *
**********************************************************
No failed systems were found in this run.
time mass zams_mass probability radius stellar_type
3.598268106262e+01 1.30592 8.75988 0.00193614 1.72498e-05 13
2.436983545180e+01 1.35842 10.9948 0.00144093 1.72498e-05 13
1.690157943855e+01 1.43124 13.7998 0.00107238 1.72498e-05 13
1.242397939734e+01 1.52416 17.3205 0.000798096 1.72498e-05 13
9.758088010684e+00 1.6691 21.7394 0.000593966 1.72498e-05 13
8.484238931453e+00 1.66351 27.2857 0.000442046 1.72498e-05 13
7.699271848690e+00 1.64817 34.247 0.000328983 1.72498e-05 13
7.568274926111e+00 1.64805 42.9844 0.000244839 1.72498e-05 13
7.569952912843e+00 1.64805 53.9508 0.000182216 1.72498e-05 13
7.571923249967e+00 1.64787 67.7151 0.00013561 1.72498e-05 13
7.612944004713e+00 1.64836 84.9909 0.000100925 1.72498e-05 13
7.620339349962e+00 1.64833 106.674 7.51114e-05 1.72498e-05 13
7.622200536471e+00 1.64833 133.89 5.59e-05 1.72498e-05 13
Binary star population
We can also set up a population that samples binary systems, by adding extra grid variables. Below is an example of a full script that runs a binary population and registers when a double compact object is formed. The logging is rather compact and should be expanded to be more useful. Also note that we run very little systems in the following example, as its just intended to show how the code works.
[14]:
def parse_function(self, output):
"""
Example parsing function
"""
# extract info from the population instance
# Get some information from the
data_dir = self.custom_options["data_dir"]
base_filename = self.custom_options["base_filename"]
# Check directory, make if necessary
os.makedirs(data_dir, exist_ok=True)
#
seperator = " "
# Create filename
outfilename = os.path.join(data_dir, base_filename)
# The header columns matching this
parameters = [
"time",
"mass_1", "zams_mass_1", "mass_2", "zams_mass_2",
"stellar_type_1", "prev_stellar_type_1", "stellar_type_2", "prev_stellar_type_2",
"metallicity", "probability"
]
# Go over the output.
for el in output_lines(output):
headerline = el.split()[0]
# CHeck the header and act accordingly
if headerline == "EXAMPLE_DCO":
values = el.split()[1:]
if not len(parameters) == len(values):
print("Number of column names isnt equal to number of columns")
raise ValueError
if not os.path.exists(outfilename):
with open(outfilename, "w") as f:
f.write(seperator.join(parameters) + "\n")
with open(outfilename, "a") as f:
f.write(seperator.join(values) + "\n")
# Create population object
example_pop = Population()
# If you want verbosity, set this before other things
example_pop.set(verbosity=0)
# Setting values can be done via .set(<parameter_name>=<value>)
example_pop.set(
# binary_c physics options
M_1=10, # bse_options
separation=0, # bse_options
orbital_period=45000000080, # bse_options
max_evolution_time=15000, # bse_options
eccentricity=0.02, # bse_options
# population_options
num_cores=2, # population_options
tmp_dir=TMP_DIR,
# Custom options: the data directory and the output filename
data_dir=os.path.join(
TMP_DIR, "example_python_population_result"
), # custom_options
base_filename="example_pop.dat", # custom_options
)
# Creating a parsing function
example_pop.set(
parse_function=parse_function, # Setting the parse function thats used in the evolve_population
)
### Custom logging
# Log the moment when the star turns into neutron
example_pop.set(
C_logging_code="""
// logger to find gravitational wave progenitors
if(stardata->star[0].stellar_type>=NS && stardata->star[1].stellar_type>=NS)
{
if (stardata->model.time < stardata->model.max_evolution_time)
{
Printf("EXAMPLE_DCO %30.12e " // 1
"%g %g %g %g " // 2-5
"%d %d %d %d " // 6-9
"%g %g\\n", // 10-11
stardata->model.time, // 1
stardata->star[0].mass, //2
stardata->common.zero_age.mass[0], //3
stardata->star[1].mass, //4
stardata->common.zero_age.mass[1], //5
stardata->star[0].stellar_type, //6
stardata->previous_stardata->star[0].stellar_type, //7
stardata->star[1].stellar_type, //8
stardata->previous_stardata->star[1].stellar_type, //9
// model stuff
stardata->common.metallicity, //10
stardata->model.probability //11
);
}
/* Kill the simulation to safe time */
stardata->model.max_evolution_time = stardata->model.time - stardata->model.dtm;
}
"""
)
# Add grid variables
resolution = {"M_1": 3, "q": 3, "per": 3}
# Mass
example_pop.add_sampling_variable(
name="lnm1",
longname="Primary mass",
valuerange=[2, 150],
samplerfunc="self.const_linear(math.log(2), math.log(150), {})".format(resolution["M_1"]),
precode="M_1=math.exp(lnm1)",
probdist="self.three_part_powerlaw(M_1, 0.1, 0.5, 1.0, 150, -1.3, -2.3, -2.3)*M_1",
dphasevol="dlnm1",
parameter_name="M_1",
condition="", # Impose a condition on this grid variable. Mostly for a check for yourself
)
# Mass ratio
example_pop.add_sampling_variable(
name="q",
longname="Mass ratio",
valuerange=["0.1/M_1", 1],
samplerfunc="self.const_linear(0.1/M_1, 1, {})".format(resolution['q']),
probdist="self.flatsections(q, [{'min': 0.1/M_1, 'max': 1.0, 'height': 1}])",
dphasevol="dq",
precode="M_2 = q * M_1",
parameter_name="M_2",
condition="", # Impose a condition on this grid variable. Mostly for a check for yourself
)
#
example_pop.add_sampling_variable(
name="log10per", # in days
longname="log10(Orbital_Period)",
valuerange=[0.15, 5.5],
samplerfunc="self.const_linear(0.15, 5.5, {})".format(resolution["per"]),
precode="""orbital_period = 10** log10per
sep = calc_sep_from_period(M_1, M_2, orbital_period)
sep_min = calc_sep_from_period(M_1, M_2, 10**0.15)
sep_max = calc_sep_from_period(M_1, M_2, 10**5.5)""",
probdist="self.sana12(M_1, M_2, sep, orbital_period, sep_min, sep_max, math.log10(10**0.15), math.log10(10**5.5), -0.55)",
parameter_name="orbital_period",
dphasevol="dlog10per",
)
# Exporting of all the settings can be done with .export_all_info()
example_pop.export_all_info()
# remove the result file if it exists
if os.path.isfile(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat")):
os.remove(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat"))
# Evolve the population
example_pop.evolve()
#
with open(os.path.join(TMP_DIR, "example_python_population_result", "example_pop.dat"), 'r') as f:
output = f.read()
print("\n")
print(output)
Grid has handled 8 stars with a total probability of 0.0211592
**********************************
* Dry run *
* Total starcount is 8 *
* Total probability is 0.0211592 *
**********************************
Signalling processes to stop
****************************************************
* Process 1 finished: *
* generator started at 2023-05-18T18:01:41.331634 *
* generator finished at 2023-05-18T18:01:42.635925 *
* total: 1.30s *
* of which 1.21s with binary_c *
* Ran 4 systems *
* with a total probability of 0.0104327 *
* This thread had 0 failing systems *
* with a total failed probability of 0 *
* Skipped a total of 0 zero-probability systems *
* *
****************************************************
****************************************************
* Process 0 finished: *
* generator started at 2023-05-18T18:01:41.321336 *
* generator finished at 2023-05-18T18:01:42.972416 *
* total: 1.65s *
* of which 1.57s with binary_c *
* Ran 4 systems *
* with a total probability of 0.0107265 *
* This thread had 0 failing systems *
* with a total failed probability of 0 *
* Skipped a total of 0 zero-probability systems *
* *
****************************************************
**********************************************************
* Population-d7df6320c8a84cbf9047b6942dd7e63b finished! *
* The total probability is 0.0211592. *
* It took a total of 2.13s to run 8 systems on 2 cores *
* = 4.27s of CPU time. *
* Maximum memory use 351.258 MB *
**********************************************************
No failed systems were found in this run.
time mass_1 zams_mass_1 mass_2 zams_mass_2 stellar_type_1 prev_stellar_type_1 stellar_type_2 prev_stellar_type_2 metallicity probability
1.406040509461e+01 1.56555 50.9713 1.65242 12.8178 13 13 13 8 0.02 0.000339963
1.817528031081e+01 1.64772 50.9713 1.41427 12.8178 13 13 13 5 0.02 0.000193036
7.588830659570e+00 1.65261 50.9713 1.65365 38.2535 13 13 13 8 0.02 0.000193036
1.199996656812e+01 1.57026 50.9713 1.64761 38.2535 13 13 13 8 0.02 0.000339963